%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reproduces Figure 1
%
% Cloyne, Dimsdale and Hürtgen (2024) 								
% "Are Tax Cuts Contractionary at the Zero Lower Bound? 
%  Evidence from a Century of Data"  		
%																			
% James Cloyne, Nicholas Dimsdale and Patrick Hürtgen, 2024
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all
clc

%% set path and solve model using dynare and Occbin

set(0,'DefaultLineLineWidth',2)
addpath('c:/dynare/5.5/matlab')

dynare Occbin_NK_baseline console

%% generate draws from prior distribution

draw2=2000;         % number of actual draws
draws=2100;         % extra number of draws in case some draws are not valid

draw_i=1;
go_flag=0;
go = true;

[beta1,beta2]=betaAB(0.5,0.2); 
[gam1,gam2]  =gammaAB(2,0.5); 
[beta5,beta6]=betaAB(0.8,0.1); 

while draw_i < draws && go

    display(['Draw ', num2str(draw_i)]);

    omegap=random('Beta',beta1,beta2);
    phipi=random('Normal',1.5,0.2); 
    while phipi<1.05
        phipi=random('Normal',1.5,0.2);
    end

    phiy=random('Normal',0.125,0.05); 
    while phiy<0
        phiy=random('Normal',0.125,0.05); 
    end

    rhot1=random('Beta',beta5,beta6);
    gamma  = random('Normal',1.5,0.25);
    adj    = random('Normal',6,1.5);
    xi     = random('gamma',gam1,gam2);

    actualdraws.omegap(draw_i)=omegap;
    actualdraws.phipi(draw_i)=phipi;
    actualdraws.phiy(draw_i)=phiy;
    actualdraws.gamma(draw_i)=gamma;
    actualdraws.xi(draw_i)=xi;
    actualdraws.adj(draw_i)=adj;
    actualdraws.rhot1(draw_i)=rhot1;

    draw_i=draw_i+1;
end

%% start priod predictive loop

tic
draw_i=1;
draws_skipped=[];
iter2=1;
iter3=1;

while iter3 < draw2

    drop_draw=0;

    omegap= actualdraws.omegap(draw_i);
    phiy  = actualdraws.phiy(draw_i);
    phipi = actualdraws.phipi(draw_i);
    xi    = actualdraws.xi(draw_i);
    gamma = actualdraws.gamma(draw_i);
    adj  = actualdraws.adj(draw_i);
    rhot1= actualdraws.rhot1(draw_i);

    set_param_value('rhot1' ,rhot1);
    set_param_value('omegap' ,omegap);
    set_param_value('phipi' ,phipi);
    set_param_value('phiy' ,phiy);
    set_param_value('adj' ,adj);
    set_param_value('gamma' ,gamma);
    set_param_value('xi' ,xi);

    options_.occbin.simul.SHOCKS=shock_mat.matrix_1;

    [oo_, out]  = occbin.solver(M_, oo_, options_);
    if out.error_flag %solution was not successful
        drop_draw=1;
    end

    options_.occbin.simul.SHOCKS=shock_mat.matrix_2;

    [oo3_, out]  = occbin.solver(M_, oo_, options_);
    if out.error_flag %solution was not successful
        drop_draw=1;
    end

    if drop_draw==0
      
        IRF.linear_Y=oo_.occbin.simul.linear(1:50+0,strmatch('y',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('y',M_.endo_names,'exact'));
        IRF.piecewise_Y=oo_.occbin.simul.piecewise(1:50+0,strmatch('y',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('y',M_.endo_names,'exact'));

        IRF.linear_INFC=oo_.occbin.simul.linear(1:50+0,strmatch('infc',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('infc',M_.endo_names,'exact'));
        IRF.piecewise_INFC=oo_.occbin.simul.piecewise(1:50+0,strmatch('infc',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('infc',M_.endo_names,'exact'));

        IRF.linear_R=oo_.occbin.simul.linear(1:50+0,strmatch('r',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('r',M_.endo_names,'exact'));
        IRF.piecewise_R=oo_.occbin.simul.piecewise(1:50+0,strmatch('r',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('r',M_.endo_names,'exact'));

        IRF.linear_TT=oo_.occbin.simul.linear(1:50+0,strmatch('tt',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('tt',M_.endo_names,'exact'));
        IRF.piecewise_TT=oo_.occbin.simul.piecewise(1:50+0,strmatch('tt',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('tt',M_.endo_names,'exact'));

        IRF.linear_RR=oo_.occbin.simul.linear(1:50+0,strmatch('rr',M_.endo_names,'exact'))-oo3_.occbin.simul.linear(1:50+0,strmatch('rr',M_.endo_names,'exact'));
        IRF.piecewise_RR=oo_.occbin.simul.piecewise(1:50+0,strmatch('rr',M_.endo_names,'exact'))-oo3_.occbin.simul.piecewise(1:50+0,strmatch('rr',M_.endo_names,'exact'));

        Simulation.y1(:,draw_i) = 100*IRF.linear_Y;
        Simulation.y2(:,draw_i) = 100*IRF.piecewise_Y;

        Simulation.r1(:,draw_i) = 400*IRF.linear_R;
        Simulation.r2(:,draw_i) = 400*IRF.piecewise_R;

        Simulation.rr1(:,draw_i) = 400*IRF.linear_RR;
        Simulation.rr2(:,draw_i) = 400*IRF.piecewise_RR;

        Simulation.invfc1(:,draw_i) = 100*IRF.linear_INFC;
        Simulation.invfc2(:,draw_i) = 100*IRF.piecewise_INFC;

        Simulation.tt1(:,draw_i) = 100*IRF.linear_TT;
        Simulation.tt2(:,draw_i) = 100*IRF.piecewise_TT;

        iter3=iter3+1;

    elseif drop_draw==1
        draws_skipped(iter2)=draw_i;
        iter2=iter2+1;
        drop_draw=0;
    end
    draw_i=draw_i + 1
end

toc

%% generate figures

start1=5;
end1=5+12;

xmin=0;
xmax=13;
stepsize=1;


h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(1);
fanChart(0:size(Simulation.y2(start1:end1,:),1)-1, Simulation.y2(start1:end1,:)); hold on;
%title('GDP at ZLB')
set(gca,'TickDir','in'); 
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('Figure1BaselineNK_GDP_ZLB');
print(filepath1, '-dpdf', '-fillpage'); print(filepath1, '-depsc') 


h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(2);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.invfc2(start1:end1,:));hold on;
%title('Inflation at the ZLB')
set(gca,'TickDir','in'); 
xlabel('Quarters')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('Figure1BaselineNK_INFL_ZLB');
print(filepath1, '-dpdf', '-fillpage'); print(filepath1, '-depsc') 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(3);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.rr2(start1:end1,:));hold on;
%title('Real rate at ZLB')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
set(gca,'TickDir','in'); 
xlabel('Quarters')
set(gca,'YLIM', [-10 80]);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('Figure1BaselineNK_RR_ZLB');
print(filepath1, '-dpdf', '-fillpage'); print(filepath1, '-depsc') 

h=gcf;
set(h,'Position',[50 50 1200 800]);
myplot=figure(4);
fanChart(0:size(Simulation.invfc1(start1:end1,:),1)-1, Simulation.tt2(start1:end1,:));hold on;
%title('Marginal tax rate away from ZLB')
plot([min(xlim()),max(xlim())],[0,0], 'k','LineWidth',0.25);
set(gca,'TickDir','in'); 
xlabel('Quarters')
set(gca,'YLIM', [-1.5 0.5]);
box off; grid off
set(gca,'XLIM', [-inf inf], 'gridlinestyle', '-.', 'gridalpha', 0.5, 'FontSize',18)
set(gca,'XLIM', [xmin xmax-1]);
set(gca,'XTick',[xmin:stepsize:xmax-1]);
myplot.PaperPositionMode = 'manual';
set(myplot,'Position',[50 50 1200 800]);
orient(myplot,'landscape')
filepath1 = fullfile('Figure1BaselineNK_TT_ZLB');
print(filepath1, '-dpdf', '-fillpage'); print(filepath1, '-depsc') 

%% compute probability of output expansion at various horizons

PROB.GDPpositive_0 =100*sum(Simulation.y1(5,:)>0)/length(Simulation.y1(5,:));
PROB.GDPpositive_1 =100*sum(Simulation.y1(6,:)>0)/length(Simulation.y1(6,:));
PROB.GDPpositive_2 =100*sum(Simulation.y1(7,:)>0)/length(Simulation.y1(7,:));
PROB.GDPpositive_3 =100*sum(Simulation.y1(8,:)>0)/length(Simulation.y1(8,:));
PROB.GDPpositive_4 =100*sum(Simulation.y1(9,:)>0)/length(Simulation.y1(9,:));
PROB.GDPpositive_5 =100*sum(Simulation.y1(10,:)>0)/length(Simulation.y1(10,:));

PROB.GDPpositive_ZLB_0 =100*sum(Simulation.y2(5,:)>0)/length(Simulation.y2(5,:));
PROB.GDPpositive_ZLB_1 =100*sum(Simulation.y2(6,:)>0)/length(Simulation.y2(6,:));
PROB.GDPpositive_ZLB_2 =100*sum(Simulation.y2(7,:)>0)/length(Simulation.y2(7,:));
PROB.GDPpositive_ZLB_3 =100*sum(Simulation.y2(8,:)>0)/length(Simulation.y2(8,:));
PROB.GDPpositive_ZLB_4 =100*sum(Simulation.y2(9,:)>0)/length(Simulation.y2(9,:));
PROB.GDPpositive_ZLB_5 =100*sum(Simulation.y2(10,:)>0)/length(Simulation.y2(10,:));

